% 各向异性扩散的有限体积法
% Gitee Repo

clc
clear

global n
global dx
global D

L = 2;
dx = 0.05;
dt = 0.0005;

[y,x] = meshgrid(-L:dx:L);
n = size(x,1);

u0 = zeros(n,n);
u1 = zeros(n,n);

u0 = exp(-2*(x.^2+y.^2));

% 各向异性扩散矩阵
% 应该是对称的
D = [0.8,0.4;0.4,0.8];

figure
clf
hold on
axis equal
h = surf(x,y,u0,'edgecolor','none');
caxis([0,1])
axis([-L,L,-L,L,0,1])
view([30,45])
drawnow

function [jx,jy] = compute_j(u0)
    global n
    global dx
    global D

    jx = zeros(n+1,n);
    jy = zeros(n,n+1);
    _jx = zeros(n+1,n);
    _jy = zeros(n,n+1);

    _jx(2:n,:) = -(u0(2:n,:) - u0(1:n-1,:)) / dx;
    _jy(:, 2:n) = -(u0(:,2:n) - u0(:,1:n-1)) / dx;

    _j = [_jx(2:n+1,:)(:)';_jy(:,2:n+1)(:)'];

    j = D*_j;

    jx(2:n+1,:) = reshape(j(1,:),n,n);
    jy(:,2:n+1) = reshape(j(2,:),n,n);

    jx(1,:)=0;
    jx(n+1,:)=0;
    jy(:,1)=0;
    jy(:,n+1)=0;
    return
end

for tick = 1:1000
    t = tick*dt;

    [jx,jy] = compute_j(u0);

    u1 = u0 - dt/dx*((jx(2:n+1,:)-jx(1:n,:))+(jy(:,2:n+1)-jy(:,1:n)));
    u0 = u1;

    if mod(tick,10) == 0
        set(h,'ZData',u1);
        set(h,'CData',u1);
        disp(sum(u1(:)))
        drawnow
        pause(0.1)
    end

end

